{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fdb5ffea",
   "metadata": {},
   "source": [
    "## Notebook 4 - Analysis of MS spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa724fbe",
   "metadata": {},
   "source": [
    "This notebook identifies products based on the database of MS2 spectra generated previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95723f60",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ../common.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f1026e25",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Determination of memory status is not supported on this \n",
      " platform, measuring for memoryleaks will never fail\n"
     ]
    }
   ],
   "source": [
    "from pyopenms import MSExperiment, MzMLFile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70a5f392",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We already prepared the reference database to analyze the main screen\n",
    "df_database = pd.read_pickle('../screening_notebooks/results/MS2_database_shifts_162_320_VB_clean.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cfd628be",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_c_glucoside(mz_peaks, i_peaks, charge, mz_precursor, adduct_type, tol=0.05, i_threshold=0.10):\n",
    "    shifts = np.array([120.042, 90.032])\n",
    "    if charge == 0:\n",
    "        if adduct_type in ['M+Glu+H','M+Glu+Na','M+Glu+NH4','M+Glu+ACN+H','M+2Glu+H','M+2Glu+Na','M+2Glu+NH4','M+2Glu+ACN+H']:\n",
    "            pass\n",
    "        elif adduct_type in ['M+Glu+2H','M+2Glu+2H']:\n",
    "            shifts = shifts / 2.\n",
    "    elif charge ==1:\n",
    "        if adduct_type in ['M+Glu+H', 'M+2Glu+H']:\n",
    "            pass\n",
    "        elif adduct_type in ['M+Glu+2H','M+Glu+Na','M+Glu+NH4','M+Glu+ACN+H','M+2Glu+2H','M+2Glu+Na','M+2Glu+NH4','M+2Glu+ACN+H']:\n",
    "            shifts = shifts / 2.\n",
    "            \n",
    "    i_peaks = i_peaks / max(i_peaks)\n",
    "    \n",
    "    for shift in shifts:\n",
    "        test = np.where(np.abs(mz_precursor - mz_peaks - shift) < tol)[0]\n",
    "        for i in test:\n",
    "            if i_peaks[i] > i_threshold:\n",
    "                return True\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac39e267",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_MS1_row(row, exp, verbose):\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"Product {row['Name']} identified with score {row['CosineScore']}\")\n",
    "\n",
    "    rt_hit = row['RT']\n",
    "    mz_hit = row['PrecursorMZ']\n",
    "\n",
    "    AUC = []\n",
    "    t_integrate = []\n",
    "    mz_ub = mz_hit*(1+ppm_integrate/1e6)\n",
    "    mz_lb = mz_hit*(1-ppm_integrate/1e6)\n",
    "    rt_lb = rt_hit - rt_semiwindow_integrate\n",
    "    rt_ub = rt_hit + rt_semiwindow_integrate\n",
    "    for spec in exp:\n",
    "        if spec.getMSLevel() == 1:\n",
    "            rt = spec.getRT()\n",
    "            if rt < rt_lb:\n",
    "                continue\n",
    "            if rt > rt_ub:\n",
    "                break\n",
    "\n",
    "\n",
    "            # Let us integrate along the m/z dimension to obtain the AUC for that time slice\n",
    "            mz_peaks, i_peaks = spec.get_peaks()\n",
    "            idx_bool = (mz_peaks > mz_lb) & (mz_peaks < mz_ub)\n",
    "            if any(idx_bool):\n",
    "                i_peaks_integrate = np.concatenate(([0.], i_peaks[idx_bool], [0.]))\n",
    "                mz_peaks_integrate = np.concatenate(([mz_lb], mz_peaks[idx_bool], [mz_ub]))\n",
    "                AUC_i = np.trapz(y=i_peaks_integrate, x=mz_peaks_integrate)\n",
    "\n",
    "                # Let us compute the baseline AUC to be substracted\n",
    "                #idx_where = np.where(idx_bool)[0]\n",
    "\n",
    "                #idx_first = idx_where[0]\n",
    "                #idx_last = idx_where[-1]\n",
    "                #i_peaks_integrate = np.array([0., i_peaks[idx_first], i_peaks[idx_last], 0.])\n",
    "                #mz_peaks_integrate = np.array([mz_lb, mz_peaks[idx_first], mz_peaks[idx_last], mz_ub])\n",
    "                #AUC_base = np.trapz(y=i_peaks_integrate, x=mz_peaks_integrate)\n",
    "                AUC_base = 0.\n",
    "\n",
    "                t_integrate.append(rt)\n",
    "                AUC.append(AUC_i-AUC_base)\n",
    "\n",
    "    # Let us integrate along the time dimension\n",
    "    t_integrate = np.concatenate(([rt_lb], t_integrate, [rt_ub]))\n",
    "    AUC = np.concatenate(([0.], AUC, [0.]))\n",
    "    AUC = np.trapz(y=AUC, x=t_integrate)\n",
    "        \n",
    "    return AUC\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32078b63",
   "metadata": {},
   "source": [
    "What we will do is to record all hits with `CosineScore` > 0.5. Then, we can plot the number of hits vs. the desired cutoff (say CosineScore = 0.85)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d16414e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "score_threshold = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe825e13",
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_mix(rt_dict={}, verbose=False):\n",
    "    \n",
    "    ### Prepare reference\n",
    "    df_mix = df_database #[df_database['Mix'] == mix_no]\n",
    "\n",
    "    ref_name = list(df_mix['Name'])\n",
    "    ref_smiles = list(df_mix['CSMILES'])\n",
    "    ref_substrates_charges = np.array(df_mix['M_charge'])\n",
    "    colnames = [a for a in df_mix.columns if 'M+' in a]\n",
    "    ref_precursors_mz = np.array(df_mix[colnames])\n",
    "    ref_fragments_mz = list(df_mix['mz'])\n",
    "    ref_fragments_i = list(df_mix['relint'])\n",
    "    ref_mona_no = list(df_mix['DB#'])\n",
    "    \n",
    "    ###\n",
    "    \n",
    "    files = glob.glob(f'../data/validation_data/mzML/*.mzML')\n",
    "    \n",
    "    df = pd.DataFrame(columns = ['File', 'Enzyme_id', 'Name', 'Rep', 'CSMILES', 'PrecursorMZ', 'Adduct', 'RT', \n",
    "                                 'CosineScore', 'C-gly', 'AUC'])\n",
    "    c = 0 \n",
    "    for filename in files:\n",
    "        c+=1\n",
    "        file = os.path.basename(filename)\n",
    "        print(f'\\n{file} (file {c} out of {len(files)})')\n",
    "\n",
    "        enzyme_id = file.split('_')[0]\n",
    "        \n",
    "        exp = MSExperiment()\n",
    "        MzMLFile().load(filename, exp)\n",
    "\n",
    "\n",
    "        ## SUBSET EXPERIMENTAL MS2 SPECTRA\n",
    "\n",
    "        list_hit_name = []\n",
    "        list_hit_mzprecursor = []\n",
    "        list_hit_rt = []\n",
    "        list_hit_score = []\n",
    "        list_hit_smiles = []\n",
    "        list_adducts = []\n",
    "        list_c_glu = []\n",
    "        list_mona_no = []\n",
    "        for spec in exp:\n",
    "            if spec.getMSLevel() == 2:\n",
    "                mz_precursor = spec.getPrecursors()[0].getMZ()\n",
    "                mz_peaks, i_peaks = spec.get_peaks()\n",
    "                rt = spec.getRT() # minutes\n",
    "                \n",
    "                test1 = ((np.abs(mz_precursor - ref_precursors_mz)/mz_precursor)*1e6) < ppm # boolean 2D matrix\n",
    "                \n",
    "                score_best = score_threshold\n",
    "                for i,j in zip(*np.where(test1)):\n",
    "                    \n",
    "                    # We also test for retention time\n",
    "                    name = ref_name[i]\n",
    "                    ref_rt = rt_dict.get(name, None)\n",
    "                    if ref_rt is None: # If not have a reference RT for a compound, assume any value is possible\n",
    "                        test2 = True\n",
    "                    else:\n",
    "                        test2 = np.abs(rt - ref_rt) < rt_semiwindow_cutoff # boolean scalar\n",
    "                    \n",
    "                    if test2:\n",
    "                        score = cosine_greedy_score(mz_peaks, i_peaks, ref_fragments_mz[i], ref_fragments_i[i])\n",
    "                        if score > score_best:\n",
    "                            score_best = score\n",
    "                            I = i\n",
    "                            J = j\n",
    "                \n",
    "                if score_best > score_threshold:\n",
    "                    list_hit_mzprecursor.append(mz_precursor)\n",
    "                    \n",
    "                    list_hit_rt.append(rt)\n",
    "                    list_hit_score.append(score_best)\n",
    "                    list_hit_smiles.append(ref_smiles[I])\n",
    "                    list_mona_no.append(ref_mona_no[I])\n",
    "                    adduct_type = colnames[J]\n",
    "                    list_adducts.append(adduct_type)\n",
    "                    name = ref_name[I]\n",
    "                    if '2Glu' in adduct_type:\n",
    "                        name += ' (2GLC)'\n",
    "                    list_hit_name.append(name)\n",
    "                    \n",
    "                    # We consider that the hit may be a C-glycoside\n",
    "                    # if it has a -120 or -90 Da neutral loss\n",
    "                    charge = ref_substrates_charges[I]\n",
    "                    is_c_glu = check_c_glucoside(mz_peaks, i_peaks, charge, mz_precursor, adduct_type)\n",
    "                    list_c_glu.append(is_c_glu)\n",
    "                    \n",
    "        \n",
    "        df_temp = pd.DataFrame({\n",
    "            'File': filename,\n",
    "            'Substrate_used': file.split('_')[1],\n",
    "            'Rep': file.split('_')[2][0],\n",
    "            'Enzyme_id': enzyme_id,\n",
    "            'Name': list_hit_name,\n",
    "            'CSMILES': list_hit_smiles,\n",
    "            'PrecursorMZ': np.round(list_hit_mzprecursor,4),\n",
    "            'Adduct': list_adducts,\n",
    "            'RT': np.round(list_hit_rt, 2),\n",
    "            'CosineScore': np.round(list_hit_score,3),\n",
    "            'C-gly': list_c_glu,\n",
    "            'MoNA_DB#': list_mona_no,\n",
    "        })\n",
    "\n",
    "        # Let us now remove duplicate products based on their name\n",
    "        # The sorting ensures we keep the highest score when dropping\n",
    "        df_temp.sort_values(by=['CosineScore'], ascending=[False], inplace=True) \n",
    "        df_temp = df_temp.drop_duplicates(subset=['Name'], keep='first') # We will keep the highest score for each Name\n",
    "        \n",
    "        print(f'{len(df_temp)} products identified')\n",
    "        \n",
    "\n",
    "        ## INTEGRATE IN THE MS1 CHANNEL (2D integration)\n",
    "        \n",
    "        list_hit_AUC = df_temp.apply(integrate_MS1_row, args=(exp, False), axis=1)\n",
    "\n",
    "        df_temp['AUC'] = np.round(list_hit_AUC,0)\n",
    "        df = pd.concat([df,df_temp], ignore_index=True)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d16a364",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Iteration 1\n",
      "\n",
      "73C5_GlycocholicAcid_1.mzML (file 1 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Osajin_3.mzML (file 2 out of 96)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/v1/xrsbvyx519g764hlvcqf01p40000gq/T/ipykernel_6102/3598349742.py:118: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.\n",
      "  df = pd.concat([df,df_temp], ignore_index=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 products identified\n",
      "\n",
      "73C5_Lobeline_3.mzML (file 3 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Steviol_3.mzML (file 4 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Boldine_3.mzML (file 5 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Capsaicin_1.mzML (file 6 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_ChlorogenicAcid_2.mzML (file 7 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_AleureticAcid_2.mzML (file 8 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_246trihydroxyacetophenone_3.mzML (file 9 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Curcumin_2.mzML (file 10 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Cortisone_2.mzML (file 11 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Esculetin_3.mzML (file 12 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Ononetin_2.mzML (file 13 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Colforsin_1.mzML (file 14 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Formononetin_3.mzML (file 15 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Formononetin_2.mzML (file 16 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Ononetin_3.mzML (file 17 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Esculetin_2.mzML (file 18 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Cortisone_3.mzML (file 19 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Curcumin_3.mzML (file 20 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Nagingin_1.mzML (file 21 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Icariin_1.mzML (file 22 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_246trihydroxyacetophenone_2.mzML (file 23 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_AleureticAcid_3.mzML (file 24 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_ChlorogenicAcid_3.mzML (file 25 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Boldine_2.mzML (file 26 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Steviol_2.mzML (file 27 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Solasodine_1.mzML (file 28 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Lobeline_2.mzML (file 29 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Osajin_2.mzML (file 30 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_18alphaglycyrrhetinicacid_2.mzML (file 31 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Galangin_3.mzML (file 32 out of 96)\n",
      "11 products identified\n",
      "\n",
      "73C5_Emodin_2.mzML (file 33 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Mevastatin_1.mzML (file 34 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Tryptophan_2.mzML (file 35 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Ornithine_1.mzML (file 36 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_10hydroxycampothecin_3.mzML (file 37 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Nagingenin_2.mzML (file 38 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Serotonin_2.mzML (file 39 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_DiffractaicAcid_3.mzML (file 40 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_DiffractaicAcid_2.mzML (file 41 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Serotonin_3.mzML (file 42 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Nagingenin_3.mzML (file 43 out of 96)\n",
      "5 products identified\n",
      "\n",
      "73C5_Emetine_1.mzML (file 44 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_10hydroxycampothecin_2.mzML (file 45 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Lagochiline_1.mzML (file 46 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Tryptophan_3.mzML (file 47 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Silibinin_1.mzML (file 48 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Galangin_2.mzML (file 49 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Emodin_3.mzML (file 50 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Quercetin_1.mzML (file 51 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_18alphaglycyrrhetinicacid_3.mzML (file 52 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Lagochiline_2.mzML (file 53 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Silibinin_2.mzML (file 54 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Mevastatin_3.mzML (file 55 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Galangin_1.mzML (file 56 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Quercetin_2.mzML (file 57 out of 96)\n",
      "11 products identified\n",
      "\n",
      "73C5_DiffractaicAcid_1.mzML (file 58 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Emetine_2.mzML (file 59 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_10hydroxycampothecin_1.mzML (file 60 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Ornithine_3.mzML (file 61 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Ornithine_2.mzML (file 62 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Emetine_3.mzML (file 63 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Nagingenin_1.mzML (file 64 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Serotonin_1.mzML (file 65 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_18alphaglycyrrhetinicacid_1.mzML (file 66 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Quercetin_3.mzML (file 67 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Emodin_1.mzML (file 68 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Mevastatin_2.mzML (file 69 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Tryptophan_1.mzML (file 70 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Silibinin_3.mzML (file 71 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_Lagochiline_3.mzML (file 72 out of 96)\n",
      "10 products identified\n",
      "\n",
      "73C5_246trihydroxyacetophenone_1.mzML (file 73 out of 96)\n",
      "6 products identified\n",
      "\n",
      "73C5_Icariin_2.mzML (file 74 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Nagingin_2.mzML (file 75 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Capsaicin_3.mzML (file 76 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Boldine_1.mzML (file 77 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Steviol_1.mzML (file 78 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Solasodine_2.mzML (file 79 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Lobeline_1.mzML (file 80 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_GlycocholicAcid_3.mzML (file 81 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Osajin_1.mzML (file 82 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Formononetin_1.mzML (file 83 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Colforsin_3.mzML (file 84 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Esculetin_1.mzML (file 85 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Ononetin_1.mzML (file 86 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_Colforsin_2.mzML (file 87 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_GlycocholicAcid_2.mzML (file 88 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Solasodine_3.mzML (file 89 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Capsaicin_2.mzML (file 90 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_ChlorogenicAcid_1.mzML (file 91 out of 96)\n",
      "7 products identified\n",
      "\n",
      "73C5_AleureticAcid_1.mzML (file 92 out of 96)\n",
      "8 products identified\n",
      "\n",
      "73C5_Icariin_3.mzML (file 93 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Nagingin_3.mzML (file 94 out of 96)\n",
      "5 products identified\n",
      "\n",
      "73C5_Curcumin_1.mzML (file 95 out of 96)\n",
      "9 products identified\n",
      "\n",
      "73C5_Cortisone_1.mzML (file 96 out of 96)\n",
      "8 products identified\n",
      "CPU times: user 8min 20s, sys: 10.7 s, total: 8min 30s\n",
      "Wall time: 8min 33s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "rt_dict = {}\n",
    "i = 1\n",
    "while True:\n",
    "    df = pd.DataFrame()\n",
    "    print(f'\\nIteration {i}')\n",
    "\n",
    "    df_tmp = process_mix(rt_dict)\n",
    "    if not df_tmp.empty:\n",
    "        df = pd.concat([df,df_tmp])\n",
    "        \n",
    "    df.to_pickle(f'./tmp/Val_results_cosine_iteration_{i}.pkl')\n",
    "    df.to_csv(f'./tmp/Val_results_cosine_iteration_{i}.csv', index=False)\n",
    "    \n",
    "    # Let us update the rt_dict for the next iteration\n",
    "    # To compute the dictionary, we get the median RT for each substrate observed across all experiments\n",
    "    # I will only use reliable hits for computing the dictionary (i.e., CosineScore>=0.85)\n",
    "    df_tmp = df[df['CosineScore']>=0.85]\n",
    "    rt_dict_new = dict(df_tmp.groupby(['Name'])['RT'].median())\n",
    "    \n",
    "    # Let us only run 1 iteration\n",
    "    if (rt_dict_new == rt_dict) | i == 1:\n",
    "        break\n",
    "        \n",
    "    rt_dict = rt_dict_new\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9700fa6f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "project_gt",
   "language": "python",
   "name": "project_gt"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
